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Abstract: The generalization of a class of low- dissipative high order filter finite difference 
methods for long time wave propagation of shock/turbulence/combustion compressible 
viscous gas dynamic flows to compressible MHD equations for structured curvilinear grids 
has been achieved. The new scheme is shown to provide a natural and efficient way for 
the minimization of the divergence of the magnetic field numerical error. Standard diver- 
gence cleaning is not required by the present filter approach. For certain MHD test cases, 
divergence free preservation of the magnetic fields has been achieved. 


1 Introduction 

An integrated approach for the control of numerical dissipation in high order finite dif- 
ference schemes in structured curvilinear grids for the compressible Euler and Navier- 
Stokes equations has been developed and verified by the authors and collaborators 
[27, 28, 18, 29, 21]. These schemes are suitable for complex multiscale compressible 
viscous flows, especially for high speed turbulence combustion and acoustics problems. 
Standard high-resolution shock-capturing schemes are too dissipative for these types of 
flow problems. For the performance of these schemes on the aforementioned flows, see 
[27, 28, 18, 29, 21, 19, 20] and references cited therein. Basically, the scheme consists of 
sixth-order or higher non-dissipative spatial difference operators as the base scheme. To 
control the amount and types of numerical dissipation, an artificial compression method 
(ACM) indicator or multiresolution wavelets are used as sensors to adaptively limit the 
amount and to aid in the selection and/or blending of the appropriate types of numeri- 
cal dissipation to be used. This adaptive control of numerical dissipation is accomplished 
by a filter step after the completion of each full time step integration of the base scheme. 
Hereafter, we refer to these schemes as the high order ACM-filter and WAV-filter methods. 

The type of base schemes used in the high order ACM-filter and WAV-filter methods 
is divergence free preserving for the magnetohydrodynamics (MHD) equations. However, 
straightforward application of the filter step to the MHD equations will not automati- 
cally preserve the divergence free magnetic field condition. With careful modification of 
the gas dynamic scheme, the filter mechanism offers several natural and efficient alter-' 
natives (without the standard divergence cleaning procedures) for minimizing the V*B 
numerical error which are not easily attainable without additional work in the standard 
high-resolution shock-capturing schemes. 'The focus of this paper is to present a filter 
approach that exhibits divergence free preservation for certain test cases. Extensive grid 
convergence comparisons with standard high-resolution shock-capturing schemes will be 
shown. 

* Proceedings of the International Conference on High Performance Scientific Computing, March 
10-14, 2003, Hanoi, Vietnam. Part of this work was performed while the second author was a 
RIACS visiting scientist at NASA Ames Research Center. 
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2 Relevance 


This paper is concerned with the compressible MHD equations, henceforth, for ease of 
reference, referred to simply as MHD equations. Throughout the paper, the term “V • B 
numerical error” refers to the “amount of non-zero value of the discretized form of V* B of 
the underlying scheme” . The following discussion pertains to schemes involving the use of 
Riemann solvers or the eigen-structure of the MHD equations. In addition, our discussion 
is restricted to the finite difference formulation for structured grids. 

An important ingredient in our method is the use of the dissipative portion of high- 
resolution shock-capturing schemes as part of the nonlinear filters. These nonlinear filters 
involve the use of approximate Riemann solvers. We will therefore first present a new form 
of high-resolution shock- capturing schemes for the conservative MHD equations using the 
non-conservative eigensystem. 

Consider the 3-D conservative and non-conservative forms of the ideal compressible 
MHD equations in Cartesian grids, 
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where the velocity vector u = (u, v,w) T , the magnetic field vector B = (B x ,B y ,B z 
is the density, and e is the total energy. The notation B 2 ~ B 2 + By -1- B 2 is used, 
pressure is related to the other variables by 


( 2 ) 

P 

The 
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$(Bl+Bl + B*)). 


For plasmas, 7 is usually equal to 5/3 (for monatomic gases). The vector on the right 
hand side of (2) is the non-conservative portion of the MHD equations [16, 17, 24]. The 
non-conservative term is proportional to (V- B). Physically, it is zero if V- B =0 initially. 
In symbolic form, the conservative and non- conservative forms can be written as 


U t +' V • F = 0, 

U t 4- V * F = 5, 

where U is the corresponding state vector, F is the conservative inviscid flux vector tensor 
and S is the non-conservative portion of the equations in (2). The non- conservative term 
S can also be written as S = Y2i=i Ni{U)U Xi , where the notation (04, £2, £3) — (x,y, z) is 
used. The curvilinear grid formulation follows the same methodology as in [25]. 
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2.1 Divergence Condition 

The V- B = 0 condition is an initial constraint for the MHD equations, and it is not part of 
the MHD differential system. This is unlike the divergence- free condition of the velocities 
for the incompressible Euler or Navier-Stokes equations which is part of the differential 
system and is needed to close the system and must be explicitly enforced. For the MHD 
equations with V - B = 0 as the initial data, all one needs is to construct schemes with the 
discretized form of V- B on the order of the truncation error, which goes to zero when the 
grid is refined. Unfortunately, straightforward extension of existing gas dynamic schemes 
to the MHD equations does not necessarily preserve the divergence-free condition. 

Presently, there are basically two camps in solving the multi-dimensional MHD equa- 
tions; namely, that which solves the conservative form, and the one which solves the 
non-conservative symmetrizable form [10, 16, 17]. For both forms of the MHD equations, 
high-resolution shock-capturing methods suffer from the need to perform extra work to 
drive the V- B numerical error down to machine zero. The popular approaches for min- 
imizing the V-B numerical error include augmenting an extra PDE to the system [2], 
using variants of the staggered approach of K.S. Yee [31, 6, 7, 4] and using a projection 
method [32]. There is a key advantage to solving the conservative equations over the non- 
conservative equations, since the conservative form guarantees correct propagation speeds 
and locations of discontinuities. The disadvantage is that the conservative form is a non- 
strictly hyperbolic system with non-convex inviscid fluxes. There exist states (e.g., triple 
umbilic points for 1-D) for which the Jacobian of the flux of the conservative form does 
not have a complete set of eigenvectors. In this paper and companion papers [22, 30], both 
the non-conservative and conservative forms of the multidimensional compressible MHD 
system are considered. 


2.2 Conservative and Non-Conservative Formulations Involving the Use of 
Approximate Riemann Solvers 


For convenience of presentation we will describe our numerical methods for the x-flux on 
a uniform grid. A more detailed discussion can be found in [22]. Let A(TJ) denote the 
Jacobian dF/dU with the understanding that the present F and S are the x-component 
of the 3-D description above. For later discussion we write the non- conservative S term in 
the x-direction as N(U)U X . 

Gallice [9] and Cargo and Gallice [1] observed that seven of the eigenvalues and eigen- 
vectors are identical for the “conservative” Jacobian matrix A and the “non-conservative” 
Jacobian matrix (A — N). The eighth eigenvector of A of the conservative system (which is 
distinct from the non-conservative system) associated with the degenerate zero eigenvalue 
can sometimes coincide with one of the other eigenvectors, thereby, making it impossible 
to define the MHD Roe’s approximate Riemann solver in the standard way. The eigenvec- 
tors of the non-conservative Jacobian ( A — N) always form a complete basis, and can be 
obtained from analytical formulas [9, 1]. A Roe type average state was developed in [9, 1]. 

We formulate our scheme together with the Gallice form of the MHD Roe’s approximate 
Riemann solver in curvilinear grids for both the conservative and non-conservative MHD 
equations. We propose to use the non-conservative form of the eigen- decomposition but 
with the degenerate eigenvalue replaced by an entropy correction [11, 26] of what was 
supposed to be the zero eigenvalue for the conservative form (e.g., a small parameter e 
that is scaled by the largest eigenvalue of A(U)). Our rationale for doing this is that 
only the eighth eigenvector of the non-conservative form is not the same as the eighth 
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eigenvector for the conservative form. The incorrect eigenvector for the conservative form 
will be multiplied by an eigenvalue which is close to zero (the eigenvalue will not be 
exactly zero when an entropy correction is used). Thus the effect of a “false” eigenvector 
will be small. By using the eighth eigenvector of the non-conservative system instead, the 
difficulty of dealing with an incomplete set of eigenvectors for the conservative system can 
be avoided. 

The conservative filter approach, and the conservative Harten-Yee, MUSCL and the 
fifth-order WENO [12] schemes used in this paper are formed by using the non-conservative 
eigen-decomposition described above in solving the conservative MHD equation set (1). 
The non-conservative filter approach, and the non-conservative Harten-Yee, MUSCL and 
WENO schemes are just the non-conservative eigen- decomposition in solving the non- 
conservative MHD equation set (2). 


3 Description of High Order Filter Methods 

Our high order ACM-filter and WAV-filter methods consist of two stages, a divergence- 
preserving base scheme stage (not involve the use of approximate Riemann solvers) and 
a filter stage (involve the use of approximate Riemann solvers). The filter stage can be 
divergence-free preserving depending on the type of filter operator being used and the 
method of applying the filter step. In order to have a good shock-capturing capability and 
improved nonlinear stability related to spurious high frequency oscillations, the blending of 
a high order nonlinear filter and a high order linear filter were proposed in our gas dynamic 
schemes. The nonlinear filter consists of the product of an ACM or wavelet sensor and 
the nonlinear dissipative portion of a high-resolution shock-capturing scheme. The high 
order linear filter is just the centered linear dissipative operator that is compatible with 
the order of the base scheme being used. 

3.1 Divergence- Free Preserving Base Scheme Step 

The first stage of the numerical method consists of a time step by a non-dissipative high 
order spatial and high order temporal base scheme operator L (e.g., a divergence-free 
preserving sixth-order central in space and fourth-order Runge-Kutta in time), 

U* = L(ET), (3) 

where U n is the numerical solution vector at time level n. When necessary, a high order lin- 
ear numerical dissipation operator can be used. For example, a divergence-free preserving 
eighth-order linear dissipation with the sixth-order centered base scheme to approximate 
F(U) X is written as 

dF 

— « DoeFj + dAx 7 {D + D^Uj, (4) 

where D QQ is the standard sixth-order accurate centered difference operator, and D + D_ 
is the standard second-order accurate centered approximation of the second .derivative. 
The small parameter d is a scaled value in the range of 0.00001 to 0.01, depending on the 
flow problem, and has the sign which gives dissipation in the forward time direction. The 
operators are modified at boundaries in a stable way [29]. 

This highly accurate base scheme is employed to numerically preserve the divergence- 
free condition of the magnetic field (to the level of round-off error) for curvilinear grids. 
When the solution is smooth, the filter step might not be needed. Thus the use of a 
high order centered difference operator as the base scheme will perfectly preserve the 
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divergence-free condition. In this case the result will be the same, whether we solve the 
conservative system (1) or non-conservative system (2). Under a shock/shear and tur- 
bulence/combustion environment, the use of a dissipative portion of the shock-capturing 
scheme as part of the filter is necessary In this case, a possible source of violation of the 
divergence-free condition can be from the filter step. 

3.2 Adaptive Numerical Dissipation Filter Step 

After the completion of a full time step of the divergence-free preserving base scheme 
stage, the second stage is to adaptively filter the solution by the product of “an ACM 
indicator or wavelet sensor” and the “nonlinear dissipative portion of a high- resolution 
shock-capturing scheme”. The final update step after the filter can be written as (assume 
1-D for ease of illustration) 

v r i = u 'i-^ H, ,w- H Ud- (s) 

The filter numerical flux vector is 

H j+ 1/2 = ^j+l/2#i+l/2- 

Here Rj+ 1/2 is the matrix of right eigenvectors of the Jacobian of the non-conservative 
MHD flux vector (Aj +x / 2 — J^j+ 1 / 2 ) evaluated at the Gallice average state Uj +1 j 2 as 

discussed in the previous subsection. The Hj+ 1 / 2 are also evaluated from the same char- 
acteristic quantities derived from these eigenvectors using the Gallice average state based 
on the U* values of (4). Due to the fact that the base scheme step is divergence free 
preserving and does not involve the use of approximate Riemann solvers, there is no dif- 
ference in solving the conservative or non-conservative system for the filter approach. To 
reduce un-necessary computations (the non-conservative portion), the non- conservative 
filter approach only solves the conservative system on the base scheme step. Thus, the 
conservative and non-conservative filter approaches differ merely by the eighth eigenvalue. 

Denote the elements of the vector H j+ x / 2 by hj + x / 2 , l = 1, 2, 8. They have the form 

fy+1/2 = ( w )J+l/2(^i+l/2)- ( 6 ) 

Here (c ^)y +1 / 2 * s a sensor to activate the shock-capturing nonlinear filter. For example, 
( w )j'+i/ 2 is designed to be zero in regions of smooth flow and near one in regions with 
discontinuities. It varies from one grid point to another and is obtained either from a 
wavelet analysis of the solution (WAV -filter scheme), or from a gradient-based detector 
(ACM- filter scheme) [27, 28, 18, 29, 21]. The blending of nonlinear filters with high order 
linear filter is discussed in [29]. 

The dissipative portion of the nonlinear filter <t> l j+1 / 2 = 9j+i/2 “^+ 1/2 is dissipative 
portion of a high order high-resolution shock- capturing scheme for the Zth-characteristic 
wave. Here g l j+1 / 2 and b l j+l / 2 are numerical fluxes of the uniformly high order high- 
resolution scheme and a high order central scheme for the Zth characteristic, respectively. 

It is noted that ^ +1 y 2 naight not be unique since there is more than one way of ob- 
taining <j> l j + y 2 ' For the forms of the <f> l j+i/ 2 use d In the numerical experiment section, see 
[27, 28, 18, 29, 21]. For example, the form of Harten and Yee and symmetric TVD schemes 
are already in the proper form in the sense that they are written in a central differencing 
portion b l j +l / 2 and a nonlinear dissipation portion N° work is required to obtain 

4 > l j+i /2 in this case. 
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With the exception of some smooth flows using the WAV-filter scheme, the filter given 
by (6), if applied to the entire MHD system (denoted by “filter all”) normally will not 
preserve the divergence free magnetic field condition. In order to minimize the numerical 
error of the divergence-free magnetic condition, the nonlinear filter step only acts on the 
gas dynamic portion of the equations (denoted by “no filter on B”). With the divergence 
free spatial base scheme and the manner that we update the solution on the filter step, 
the divergence free property should be preserved by the “no filter on B” option. There are 
additional variants of the filter approach that from a theoretical standpoint, are divergence 
free. See [30] for more details. 


4 2-D Compressible MHD Numerical Examples 


For illustrative purposes, numerical experiments using sixth-order central spatial dis- 
cretization as the base scheme is chosen for the ACM-filter and WAV-filter schemes. The 
sixth-order base scheme together with the nonlinear /linear filter with wavelet sensor will 
be denoted WAV66. When a more conventional gradient based sensor ACM is used, the 
scheme is denoted ACM66. If high order linear numerical dissipation is also used in the 
base scheme, the methods will be denoted WAV66+AD8 and ACM66+AD8 respectively. 
The strength of the eighth-order dissipation will be denoted by a tunable coefficient d, 
as in (4). In all of the filter scheme computations, the nonlinear dissipative portion of 
Harten-Yee is used as part of the nonlinear filter term. For all test cases, the entropy fix 
parameter is 0.25 for the ACM and WAV-filter schemes. 

The fifth-order weighted ENO scheme [12] (WEN05), and second-order Harten-Yee 
and MUSCL schemes are used for comparison. Classical fourth-order Runge-Kutta time 
stepping is used for all sixth-order schemes, as well as for the WEN05 scheme. The second- 
order Harten-Yee and MUSCL are integrated in time by the second-order TVD Runge- 
Kutta method. Except for WEN05, the minmod limiter, the van Leer version of the van 
Albada limiter and the Colella- Woodward limiter are considered. 

The V- B numerical error is obtained by approximating the spatial derivatives by sixth- 
order centered differences for WAV66, ACM66 and WEN05, whereas the corresponding 
V- B numerical error is obtained by second-order centered differences for the second-order 
TVD schemes (MUSCL and Harten-Yee). The L 2 -norm of V B of a particular scheme 
is computed by taking the square root of the sum over the square of all three spatial 
directions of the discretized form of V- B at all grid points. 

4.1 MHD Kelvin-Helmholtz Instabilities (7 = 1.4, Periodic BC) 

The magnetohydro dynamic Kelvin-Helmholtz instabilities have been studied by many pre- 
vious investigators [2, 13, 8]. We have used the set up in [2] which is shown in Fig. 1. Snap- 
shots of the time evolution of the ^-velocity is also shown in Fig. 1 by CEN66-f AD8 (sixth- 
order central with an eighth-order linear dissipative added to the base scheme (d ~ 0.001). 
The solution is obtained without the filter step. At stopping time T — 0.5, the problem 
is smooth enough that it can be solved by the base scheme alone. Density contours at 
time T = 0.5 with 30 equidistant contour levels between 0.4 and 1.2 are used. Five levels 
of grid refinement are considered, namely, 51 x 101, 101 x 201, 201 x 401, 401 x 801 and 
801 x 1601. Grids of increasing refinements by the eighth-order central difference with a 
tenth-order linear dissipation added (CEN884-AD10, d = 0.001) are used as the reference 
solution. Computations using d = 0 (CEN88) are not stable for the five grids. 
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MHD Kelvin-Helmholtz Instability (k - 1 .*) 
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7 = 0.32426 T = 0.39551 T = 0.49229 



0 0.5 1 0 0.5 1 O 0.5 1 
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Time discretization: 2nd-order Runge Kutta - 2nd-order spatial scheme 

4th-order Runge-Kutta - 6th or higher-order spatial scheme 


Fig. 1. Problem setup and time evolution of the Kelvin-Helmholtz problem, x- velocity contours 
bv CEN66+AD8 on 101x201 grid points. 


MHD Kelvin-Helmholtz Instability (T=0.5) 
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Fig. 2. Density (left) and V- B (middle) contours at T — 0.5, and L 2 -norm of V- B as a function 
of time (right) by MUSCL (top row) and WEN05 (bottom row). 
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MHD Kelvin-Helmholtz Instability (T=0.5) 

(Divergence Free ACM66 & WAV66, 201 x 401) 
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Fig. 3. Density (left) and V- B (middle) contours at T — 0.5, and L 2 -n orm of V B as a function 
of time (right) by ACM66 (top row) and WAV66 (bottom row). 


At time 0.5 the problem is smooth enough and there is no need for the more CPU 
intensive shock-capturing schemes. However, as the flow evolves at a later time, shock- 
capturing methods are required. Here, the purpose is to examine the V* B numerical error 
when the flow is still smooth using shock-capturing methods. Figure 2 shows the density 
(left) and VB (middle) contours at T = 0.5, and L 2 -norm of V-B as a function of 
time (right) by MUSCL (top row) and WEN05 (bottom row). The same computations 
by ACM66 and WAV66 are shown in Fig. 3. The V-B contours with 30 equidistant 
contour levels between —150 and 150 are used. The CPU time used was considerably larger 
(around a factor 2.5) for the WEN05 scheme. MUSCL, Harten-Yee and WEN05 exhibit 
small oscillations at the outer edges of the vortices as the grid is refined. It is possible to 
decrease these oscillations by increasing the multidimensional entropy fix parameters of 
the Harten-Yee scheme [26]. 

Density contours using ACM66, ACM66+AD8, WAV66 and WAV66+AD8 for all lim- 
iters exhibit an accuracy similar to CEN884-AD8 (figures not shown). There is no gain 
in solving the conservative over the non-conservative system for these two filter schemes. 
However, their V-B numerical errors are very different when using the “no filter on B” op- 
tion verses the “filter all” option. They are also very different from the standard MUSCL, 
Harten-Yee and WEN05 schemes. For the no filter on B equations option, divergence free 
preservation is achieved by the ACM66 and WAV66. The three standard shock-capturing 
methods exhibit similar V- B numerical errors. 

4.2 A 2-D Compressible MHD Riemann Problem (7 = 5/3, Dirichlet BC) 

We examine the same 2-D Riemann problem as in [2], It consists of four constant states 
at time zero, as shown in Fig. 4. Grid convergence studies solving conservative (top) and 
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2 -D MHD Riemann Problem (r = 5/3) 


i.c. 
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Fig. 4. Schematic of the initial data for the 2-D Riemann problem. 


2 -D MHD Riemann Problem 


Con. vs. Noncon. 
WEN05, T=0.2 


WEN05 

201x201 401x401 801x801 



Fig. 5. Grid refinement by WEN05. Density contours solving the conservative (upper row) and 
non-conservative (lower row) form of the equations. 
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non-conservative (bottom) system by WEN05 are shown in Fig. 5 for density contours at 
T — 0.2 with 40 contours equally spaced between 0.75 and 2.L 



The accuracy in a solution of a Riemann problem away from discontinuities is difficult 
to improve by increasing the order of the scheme. A large part of the solution is constant, 
and the structure that develops is affected by low order errors from the discontinuity in the 
initial data. Since all five methods can capture shocks within 2-4 grid cells, their density 
contours look very similar even-though the V-B contours or the L 2 norm of the V-B 
numerical errors are all very different. 

The effect on V ■ B when switching from a non-conservative system to a conservative 
system is less significant for the Harten-Yee and WEN05 than for MUSCL. V-B con- 
tours* for the three methods, MUSCL, Harten-Yee, and WENO-5 are displayed in Fig. 6. 
The V-B contours use 30 equidistant contour levels between -3.7 and 3.7. The ACM66, 
ACM66+AD8 WAV66 and WAV66+AD8 methods all exhibit divergence free preservation 
when no nonlinear filter is applied on the B equations. Figures 7 show a comparison. 
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2 -D MHD Riemann Problem 

(Filter All) vs. (No Filter on B) 

ACM66+AD8 


\ Filter All] j No Filter on B 

|L2 - norm vs. T ime! 



0 0.1 0.2 0 0.1 0,2 


Div(B), T=0.2 

1 

0,5 
>* 0 
-0.5 

-1 0 1 ~-1 0 1 
X X 

Fig. 7. L 2 norm of V-B vs. time and V-B contours at T — 0.2 by ACM66+AD8 when the non- 
linear filter is not applied on the magnetic field (left) and when it is applied to all components 
(right). 201 x 201, 401 x 401, and 801 x 801 grid points. 

4.3 Compressible MHD Orszag-Tang Vortex (7 — 5/3, Periodic BC) 

The 2-D Compressible MHD Orszag-Tang vortex problem [14, 3, 15, 4, 5] consists of 
periodic boundary conditions with smooth initial data is shown in Fig. 8. Density contours 
by WAV66+AD8 at T — 3.14 using “filter all” and “no filter on B” are shown in the same 
figure. The density contours are almost identical. 

The initial sine waves break into discontinuities at a later time with complicated flow 
interactions. The computation stops at time T= 3.14 (w 7 r), when discontinuities have 
formed and interacted. The solution has both complicated structure and discontinuities. It 
is a problem well suited for demonstrating our approach with highly accurate methods for 
solutions with discontinuities. Density contours with 30 equally spaced contours between 
0.9 and 6.1 are used for illustration. Again, the same five levels of grid refinement study 
were performed on all five methods. 

Figures 9 and 10 show the comparison of WEN05 (solving both systems) with 
WAV66+AD8. Divergence free preservation is achieved by WAV66+AD8 using the “no 
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Compressible Orszag-Tang Vortex (y = 5 / 3 ) 



Fig. 8. Schematic, problem setup and density contours by WAV664-AD8 for the Orszag-Tang 
problem using a 801 x 801 grid at time T — 3.14. 

Compressible Orszag-Tang Vortex 

Con. vs. Noncon. 

WEN05 


801x801 

Con. Noncon. 

L 2 - norm vs. Time 



X X 

Fig. 9. L 2 norm of V- B in time (top) and V- B contours (bottom) by WEN05 for the conservative 
(left) and the non-conservative (right) systems. 
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Fig. 10. WAV66+AD8 L 2 norm of V B in time (top row), V-B contours at T = 3.14 (bottom 
row). No non-linear filter on B (left) and non-linear filter on all components (right). 

filter on B” option. ACM664-AD8 exhibits a similar behavior as WAV66+AD8 with the 
exception that divergence free is also possible for the “filter all” option for WAV66+AD8 
for T < 0.7, whereas the ACM66+AD8 losses its divergence free preservation at a much 
earlier time. The behavior of WAV66 (d ~ 0) and ACM66 ( d = 0) is similar. 

The resolution of the global structure of the density contours is well captured by all five 
methods. However, small fine structures were captured by the ACM-filter and WAV-filter 
schemes on a 101 x 101 grid, and not by MUSCL, Harten-Yee and WEN05. 

4.4 A Planar Shock Interacting with a Magnetic Cloud (7 = 5/3, Supersonic 
Inflow &; Open Boundaries) 

The fourth test problem is a planar shock interacting with a magnetic cloud studied in 
[4, 5]. This is a more challenging problem to simulate. The same initial configuration as in 
[23] is considered here. The problem setup and schematic of the initial condition is shown 
in Fig. 11. Figure 11 also shows the density contours by WEN05 at T ~ 0.06. Figures 12 
and 13 show the same comparison between WEN05 and WAV66+AD8 as the previous 
problem. The behavior of ACM66-f*AD8, WAV66-fAD8, WAV66 ( d = 0) and ACM66 
(d = 0) is similar to the Orszag-Tang problem. The same conclusion can be drawn for 
evolution time up to T = 0.04. 

This problem is very stiff and very small CFL is required. For the finer grid, in order 
to obtain a stable solution by WEN05, the CPU time is more than an order of magnitude 
greater than for the Harten-Yee and MUSCL schemes, and many times more CPU than 
the ACM and WAV-filter scheme. This is partially due to a lower stability limit of WEN05 
than the rest of the schemes. For all four test cases, MUSCL and Harten-Yee require similar 
CPU. The ACM and WAV-filter schemes require slightly more CPU than the Harten-Yee 
and MUSCL schemes. For almost all problems, WEN05 requires more CPU than ACM 
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Fig. 11. Problem setup and schematic of the initial data of the MHD shock/cloud interaction 
problem. Density contours by WEN05 solving the conservative system using a 801 x 801 grid. 


and WAV-filter schemes. In addition, for all test cases and all five methods, the V*B 
contour numerical errors (at their corresponding stopping times) increase as the grid is 
refined. For a detailed comparison and the performance of all five schemes of all the test 
cases, see [30]. 


5 Concluding Remarks 

A natural and efficient high order filter approach in the sense of not needing traditional 
divergence cleaning for the minimization of the V- B numerical error was proposed and val- 
idated using four 2-D compressible MHD test cases. Five levels of grid refinement on four 
different flow types were compared with three standard high-resolution shock-capturing 
schemes, namely, a second-order MUSCL and Harten-Yee upwind TVD schemes, and the 
fifth-order WENO scheme. The role that the proper treatment of the corresponding nu- 
merical boundary conditions can play on the effect of reducing the V • B numerical error 
was studied in [30]. Among the four test cases, with the exception of the MHD shock/cloud 
interaction problem, we can safely conclude that divergence free high order filter schemes 
for the compressible MHD equations are possible without the need of standard divergence 
cleaning. These schemes are applicable to a wide variety of flow physics problems. Appli- 
cation of these schemes to viscous MHD flows with resistivity and multiscale structure is 
forthcoming. 
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